C*************************************************************
C*                           CAS2D                           *
C*-----------------------------------------------------------*
C*  CALCULATION OF TRANSONIC INVISCID POTENTIAL FLOW THROUGH *
C* 2-D STATIONARY AIRFOIL CASCADES USING FINITE AREA METHOD, *
C*   ROTATED TYPE DEPENDENT FINITE DIFFERENCING,ARTIFICIAL   *
C*        VISCOSITY AND SUCCESSIVE LINE OVER-RELAXATION      *
C*                          -*-                              *
C*                     WINTER 1979.                          *
C*               DR. DJORDJE S. DULIKRAVICH                  *
C*    NATIONAL RESEARCH COUNCIL - NASA RESEARCH ASSOCIATE    *
C*************************************************************
C*      NASA LEWIS RESEARCH CENTER,CLEVELAND,OHIO 44135      *
C* COMPUT.FLUID MECH.BRANCH,5-9,PHONE:(216)433-4000,EXT.6859 *
C*-----------------------------------------------------------*
C*            CORNELL UNIVERSITY,ITHACA,N.Y. 14853           *
C* MECH.&AERO.ENG.DEP.,UPSON HALL,RM 216,PHONE:(607)256-3372 *
C*************************************************************
C------ PNTS ARE NUMBERED FROM LOWER T.E. TOWARDS THE L.E.
C*************************************************************
C*                           CAS2D                           *
C*-----------------------------------------------------------*
C*  CALCULATION OF TRANSONIC INVISCID POTENTIAL FLOW THROUGH *
C* 2-D STATIONARY AIRFOIL CASCADES USING FINITE AREA METHOD, *
C*   ROTATED TYPE DEPENDENT FINITE DIFFERENCING,ARTIFICIAL   *
C*        VISCOSITY AND SUCCESSIVE LINE OVER-RELAXATION      *
C*                          -*-                              *
C*                     WINTER 1979.                          *
C*               DR. DJORDJE S. DULIKRAVICH                  *
C*    NATIONAL RESEARCH COUNCIL - NASA RESEARCH ASSOCIATE    *
C*************************************************************
C*      NASA LEWIS RESEARCH CENTER,CLEVELAND,OHIO 44135      *
C* COMPUT.FLUID MECH.BRANCH,5-9,PHONE:(216)433-4000,EXT.6859 *
C*-----------------------------------------------------------*
C*            CORNELL UNIVERSITY,ITHACA,N.Y. 14853           *
C* MECH.&AERO.ENG.DEP.,UPSON HALL,RM 216,PHONE:(607)256-3372 *
C*************************************************************
C------ PNTS ARE NUMBERED FROM LOWER T.E. TOWARDS THE L.E.
C------ WITH THE PRESENT COMMON REGION, THE MAX. NUMBER OF
C---    INPUT POINTS IS 129 AND THE MAX. NUMBER OF MESH
C---    CELLS IS (192X48).
C------ NON-SYMMETRIC AIRFOIL : NUMBER OF INPUT PNTS MUST BE ODD.
C---    (MEANING THAT ALL THE POINTS MUST BE PROVIDED IN INPUT,
C---    COUNTING T.E. TWICE)
C------ SYMMETRIC AIRFOIL : NUMBER OF INPUT PNTS MUST BE EVEN.
C---    (MEANING THAT ONLY THE LOWER SURFACE POINTS SHOULD BE GIVEN)
C*-----------------------------------------------------------*
C---    XCELL        =NO. OF MESH CELLS ON THE AIRFOIL (FIRST GRID)
C---    YCELL        =NO. OF MESH CELLS OFF THE AIRFOIL (FIRST GRID)
C---    PMESH        =TOTAL NUMBER OF GRIDS USED (MAX. FOUR)
C---    TWIST        =STAGGER ANGLE (IN DEGREES)
C---    ALPHA1,ALPHA2=FREE STREAM ANGLES AT INFINITIES (IN DEGREES)
C---    PITCH        =GAP-TO-CHORD RATIO
C---    RO1,RO2      =RADII OF THE LEADING AND TRAILING EDGES
C---    FMACH        =MACH NO. AT UPSTREAM INFINITY
C---    RLX          =RELAXATION FACTOR ON THE FIRST GRID
C---    AR           =STREAMTUBE INLET AREA OVER EXIT AREA RATIO
C---    TITR1        =NUMBER OF ITERATION SWEEPS ON THE FIRST GRID
C---    TITR2        =NUMBER OF ITERATION SWEEPS ON THE SECOND GRID
C---    TITR3        =NUMBER OF ITERATION SWEEPS ON THE THIRD GRID
C---    POINTS       =NO. OF PNTS WHOSE COORD. ARE GIVEN IN INPUT
C---    GAMMA        =RATIO OF THE SPECIFIC HEATS OF THE GAS
C---    CONVER       =RELATIVE CIRCULATION CONVERGENCE RATE CRITERION
C*-----------------------------------------------------------*
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DIMENSION TITLE(20)
      DOUBLE PRECISION DSIN,DCOS,DLOG,DATAN,DARSIN,RO1,RO2,GEND,        -
     2 WV,DEXP,PITCH,DSQRT,TWIST,EM,C1C2,XOR,YOR,SLOPEX,COEF,G22,       -
     3 RLX,ALM,XSI,TOM,DABS,PP,QQ,RO,CIRC,FO,CFST,EM2,CIRCO,GORE
      DOUBLE PRECISION ALPHA1,ALPHA2,PI,RADIAN,AP,AQ,FMACH,AR,A,        -
     2 XCELL,YCELL,SUB,HR,VR,POTCH,SIS,DRD,CHORD,PMESH,AC,BC,DC,        -
     3 RUK,TW,TRI,SOS,EMM,EMP,GAMMA,CONVER,HALF,EXSIN,EXCOS,WAKE,B,     -
     4 TITR1,TITR2,TITR3,POINTS,XLE,YLE,XTE,YTE,XLEAD,YLEAD,CHORDM
      RADIAN=57.2957795130823D+00
      PI=3.14159265358979D+00
1     FORMAT(20A4)
      WRITE(6,2)
2     FORMAT(//,35X,43HCALCULATION OF TRANSONIC INVISCID POTENTIAL)
      WRITE(6,3)
3     FORMAT(35X,42HFLOW THROUGH 2-D STATIONARY CASCADE(CAS2D))
      WRITE(6,4)
4     FORMAT(35X,44HPROGRAM DEVELOPED BY: DJORDJE S. DULIKRAVICH)
5     FORMAT(2X,7X,D12.5,7X,D12.5,7X,D12.5)
C------ INPUT PARAMETERS
      READ(5,1) (TITLE(I),I=1,20)
      WRITE(6,7) (TITLE(I),I=1,20)
7     FORMAT(//,35X,20A4,/)
      READ(5,5) XCELL,YCELL,PMESH
      READ(5,5) ALPHA1,TWIST,ALPHA2
      READ(5,5) PITCH,RO1,RO2
      READ(5,5) FMACH,RLX,AR
      READ(5,5) TITR1,TITR2,TITR3
      READ(5,5) POINTS,GAMMA,CONVER
      WRITE(6,8)
8     FORMAT(10X,5HXCELL,15X,5HYCELL,15X,5HFMACH,15X,5HGAMMA,           -
     1      15X,6HCONVER)
      WRITE(6,9) XCELL,YCELL,FMACH,GAMMA,CONVER
9     FORMAT(3X,5D20.8)
      WRITE(6,10)
10    FORMAT(/,10X,6HITRMAX,15X,3HRLX,17X,3HRO1,17X,3HRO2,17X,5HPMESH)
      WRITE(6,9) TITR1,RLX,RO1,RO2,PMESH
C------ READING X&Y COORD. FROM UNIT NO.5
      MAXP=POINTS+0.0001
      READ(5,11) (HR(I),VR(I),I=1,MAXP)
11    FORMAT(2F10.6)
      MPHF=(MAXP+1)/2
      P1=POINTS+0.1
      I1=P1/2
      P2=2*I1+0.2
      IF(P1 .GT. P2) GO TO 14
      MPHF=MAXP
      MAXP=2*MPHF-1
      DO 12 II=1,MPHF
      I=MAXP+1-II
      HR(I)=HR(II)
      VR(I)=-VR(II)
12    CONTINUE
14    CONTINUE
      CHORDM=0.0D+00
      XTE=HR(1)
      YTE=VR(1)
      DO 16 I=5,MAXP
      DX=HR(I)-XTE
      DY=VR(I)-YTE
      CHORD=DX*DX+DY*DY
      IF(CHORD .LT. CHORDM) GO TO 18
      CHORDM=CHORD
16    CONTINUE
18    CONTINUE
      ILE=I-1
      CHORD=DSQRT(CHORDM)
      XLE=HR(ILE)
      YLE=VR(ILE)
      XOR=XTE-XLE
      YOR=YTE-YLE
      TW=YOR/XOR
      DRD=TWIST
      IF(DABS(TW) .LT. 0.001) GO TO 21
      TW=DATAN(TW)
      PP=DSIN(TW)
      QQ=DCOS(TW)
      TWIST=TWIST/RADIAN
      TWA=DSIN(TWIST)
      TWB=DCOS(TWIST)
      NTW=1
      IF(TWIST .LT. 0.0D+00) NTW=-1
      TWIST=TWIST-((NTW-1)/2)*PI
      DO 20 I=1,MAXP
      S(I)=HR(I)-XLE
      F(I)=VR(I)-YLE
      HR(I)=QQ*S(I)+PP*F(I)
      VR(I)=QQ*F(I)-PP*S(I)
      HR(I)=HR(I)/CHORD
      VR(I)=VR(I)/CHORD
      S(I)=TWB*HR(I)+TWA*VR(I)
      F(I)=TWB*VR(I)-TWA*HR(I)
20    CONTINUE
      PP=DSIN(TWIST)
      QQ=DCOS(TWIST)
      GO TO 25
21    CONTINUE
C----- AIRFOIL GIVEN AT ZERO TWIST ANGLE
      TWIST=TWIST/RADIAN
      TWA=DSIN(TWIST)
      TWB=DCOS(TWIST)
      NTW=1
      IF(TWIST .LT. (-1.0D-05)) NTW=-1
      TWIST=TWIST-((NTW-1)/2)*PI
      PP=DSIN(TWIST)
      QQ=DCOS(TWIST)
      DO 22 I=5,MAXP
      IF(HR(I) .GT. HR(I-1)) GO TO 23
22    CONTINUE
23    XLE=HR(I-1)
      YLE=VR(I-1)
      DO 24 I=1,MAXP
      HR(I)=(HR(I)-XLE)/CHORD
      VR(I)=(VR(I)-YLE)/CHORD
      S(I)=TWB*HR(I)+TWA*VR(I)
      F(I)=TWB*VR(I)-TWA*HR(I)
24    CONTINUE
25    CONTINUE
      TWIST=TWIST+((NTW-1)/2)*PI
      TW=TWIST*RADIAN
      WRITE(6,26)
26    FORMAT(/,10X,6HALPHA1,15X,5HTWIST,15X,6HALPHA2,15X,2HAR,          -
     1       17X,5HPITCH)
      WRITE(6,9) ALPHA1,TW,ALPHA2,AR,PITCH
      WRITE(6,27)
27    FORMAT(/,21X,21HAIRFOIL AT ZERO ANGLE,18X,                        -
     2 30HAIRFOIL AT CORRECT INCLINATION)
      WRITE(6,28) (HR(I),VR(I),S(I),F(I),I,I=1,MAXP)
28    FORMAT(5X,4F20.5,I10)
      MAXX=XCELL+3.001
      MAXXM=MAXX-1
      AC=DSIN(TWIST)
      BC=DCOS(TWIST)
      ALPHA1=ALPHA1/RADIAN
      AP=DSIN(ALPHA1)
      AQ=DCOS(ALPHA1)
      ALPHA2=ALPHA2/RADIAN
      EXSIN=DSIN(ALPHA2)
      EXCOS=DCOS(ALPHA2)
      IMID=(MAXX+1)/2
      MAXY=YCELL+2.001
      MAXYM=MAXY-1
      MAXYP=MAXY+1
      NMESH=PMESH*1.0001
      ITRMAX=TITR1+0.001
      TITR=TITR1
C------ ITERATIVE DETERMINATION OF THE DOWNSTREAM INFINITY B.C.
      TRI=GAMMA-1.0D+00
      RUK=1.0D+00/TRI
      TOM=0.9
      SOS=0.5D+00*TRI*FMACH*FMACH
      XSI=EXCOS/(AQ*AR)
      IT=1
30    CONTINUE
      YOR=1.0D+00+SOS*(1.0D+00-TOM*TOM)
      RO=YOR**RUK
      SIS=TOM*XSI-1.0D+00/RO
      DC=XSI*(YOR-TOM*TOM*FMACH*FMACH)
      TOM=TOM-YOR*SIS/DC
      IF(IT .GT. 200) GO TO 800
      IT=IT+1
      IF(DABS(SIS) .GT. 1.0D-12) GO TO 30
      RO=TOM*XSI
      RO=1.0D+00/RO
      CFST=PITCH*(AP-TOM*EXSIN/AR)
      XOR=(1.0D+00/(FMACH*FMACH))+0.5D+00*TRI*(1.0D+00-TOM*TOM)
      XOR=TOM/DSQRT(XOR)
      YOR=0.5D+00*RO*(3.0D+00-TOM*TOM)
      WRITE(6,32)
32    FORMAT(/,40X,30HDOWNSTREAM INFINITY PARAMETERS)
      WRITE(6,34)
34    FORMAT(20X,8HVELOCITY,7X,7HDENSITY,7X,8HMACH NO.,7X,              -
     1       11HCIRCULATION,5X,8HPRESSURE)
      WRITE(6,36) TOM,RO,XOR,CFST,YOR
36    FORMAT(14X,5E15.6)
      POTCH=0.5D+00*PITCH
      SLOPEX=PI/POTCH
      NHNOW=1
      C1C2=(1.0D+00-(RO1+RO2)*0.5D+00)*SLOPEX
      ALM=0.5D+00*C1C2+RO1*0.5D+00*SLOPEX+PI*PP
C------ ITERATIVE DETERMINATION OF THE SLIT LENGTH IN Z-PLANE
      EM=1./DEXP(PITCH)+0.1
      N=1
38    CONTINUE
      IF(DABS(EM) .GT. 0.999) EM=0.95
      EM2=EM*EM
      EMP=1.0D+00+EM2
      EMM=1.0D+00-EM2
      HALF=2.0D+00*EM*QQ/EMM
      XSI=1.0D+00+HALF*HALF
      XOR=HALF+DSQRT(XSI)
      XOR=QQ*DLOG(XOR)
      YOR=2.0D+00*EM*PP/EMP
      YOR=PP*DARSIN(YOR)
      SIS=RO1*0.5D+00*SLOPEX-ALM
      A=-PI*PP-2.0D+00*(XOR+YOR)-SIS
      XOR=EMP/DSQRT(XSI)
      XOR=XOR/(EMM*EMM)
      YOR=EMP*EMP-4.0D+00*EM2*PP*PP
      YOR=EMM/(EMP*DSQRT(YOR))
      B=-4.0D+00*(QQ*QQ*XOR+PP*PP*YOR)
      YOR=EM-A/B
      XOR=YOR-EM
      EM=YOR
      N=N+1
      IF(N .GT. 200) STOP
      IF(DABS(XOR) .GT. 1.D-15) GO TO 38
C------ ANGLE OF THE TRAILING EDGE IN THE CIRCLE PLANE
      EM2=EM*EM
      FO=(1.0D+00-EM2)/(1.0D+00+EM2)
      FO=FO*PP/QQ
      FO=DATAN(FO)+((NTW+1)/2)*PI
C------ STRETCHING,SHIFTING AND SHEARING OF THE COORDINATES
      B=PI*QQ
      DO 40 I=1,MAXP
      HR(I)=HR(I)*SLOPEX-ALM
      VR(I)=VR(I)*SLOPEX+B
40    CONTINUE
      VR(MAXP)=0.5D+00*(VR(1)+VR(MAXP))
      VR(1)=VR(MAXP)
      HR(MAXP)=0.5D+00*(HR(1)+HR(MAXP))+1.0D-10
      HR(1)=HR(MAXP)
      CALL CONMAP(MAXP,EM,FO,PP,QQ,NTW)
      DRD=PI-0.5D+00*(S(MAXP)-S(1))
      S(1)=S(1)-DRD*NTW
      S(MAXP)=S(MAXP)+DRD*NTW
      DRD=S(MAXP)-2.0D+00*PI
C------ PERIODIC EXTENTION OF THE MAPPED REGION
      DO 42 I=2,MPHF
      II=MAXP+I-1
      S(II)=S(I)+2.0D+00*PI
      F(II)=F(I)
42    CONTINUE
      MAXP=MAXP+MPHF-1
      XOR=2.0D+00*PI+DRD
      DO 44 I=1,MAXP
      IF(S(I) .GT. XOR) GO TO 46
44    CONTINUE
46    CONTINUE
      MAXP=I-1
      CALL SPLIF(1,MAXP,3,0.,3,0.,0,0.)
C------ INTERPOLATED 'VE'-COORD. OF THE AIRFOIL SURFACE
      MX=NMESH-1
      YOR=2.0D+00**MX
      MX=XCELL*YOR+2.0001
      COEF=4.0D+00*PI/(XCELL*YOR)
      GORE=0.02-0.015/PITCH-0.5*RO1
      DO 55 I=2,MX
      YOR=COEF*(I-2)
      SI(I)=DRD+0.5D+00*YOR+GORE*DSIN(YOR)
55    CONTINUE
      CALL INTPL(2,MX,1,MAXP,0)
C------ SAVE INTERPOLATED VALUES
      DO 60 I=2,MX
      SRF(I)=FI(I)
60    CONTINUE
      SRF(2)=0.5D+00*(SRF(2)+SRF(MX))
      SRF(MX)=SRF(2)
      NON=NMESH-1
      NON=2.0D+00**NON+0.0001
      I=2
      DO 62 II=2,MX,NON
      FI(I)=SRF(II)
      I=I+1
62    CONTINUE
C------ INITIAL GUESS FOR THE REDUCED POTENTIAL FIELD
      DO 80 J=1,MAXYP
      DO 80 I=1,MAXX
      G(I,J)=0.0
80    CONTINUE
      CIRC=0.0001*CFST
      CIRCO=0.0001
      IF(NMESH .GT. 1) CONVER=0.5*CONVER/(PMESH-1.)
C------ AFTER EACH MESH REFINEMENT START FROM HERE AGAIN
100   CONTINUE
      CALL REMAP(GORE,ALM,PP,QQ,XCELL,YCELL,EM,DRD,NTW,AC,BC,RO2)
      CALL XYINF
      G(IRHS,2)=0.0D+00
      G(IMID,2)=0.0D+00
      G(ILHS,2)=0.0D+00
      GEND=(TOM*EXCOS-AQ)*(X(2,2)-X(2,3))
      MIMA=9/NHNOW
      CMIN=0.4*CFST
      MXX=MAXXM-2
      MYY=MAXY-2
      WRITE(6,105) MXX,MYY
105   FORMAT(//,32X,44HRESULTS OF THE ITERATIVE PROCEDURE ON MESH =,    -
     1       I3,1X,1HX,I3,/)
      WRITE(6,110)
110   FORMAT( 10X,4HITER,4X,2HIR,3X,2HJR,2X,12HMAX. RESIDUE,4X          -
     1,2HIG,3X,2HJG,3X,12HMAX.CORRECT.,3X,11HCIRCULATION,2X,            -
     212HRELAX. COEF.,1X,4HISTG,2X,4HNSUP)
C------ ITERATION PROCEDURE FOR POTENTIAL FIELD STARTS HERE
      DO 600 ITER=1,ITRMAX
      SUB=RLX
      IF(DABS(CIRC) .LT. CMIN .OR. ITER .LE. MIMA) SUB=1.
      CIRCO=CIRC
      XOR=((ITER-1)*1.0D+00)/TITR
      CIRC=CIRC+(CFST-CIRC)*XOR
      YOR=0.
      IF(DABS(CFST) .GT. 1.0D-06) YOR=(CIRC-CFST)/CFST
      IF(DABS(YOR) .LT. 0.05) CIRC=0.5D+00*(CFST+CIRC)
      G22=GEND+0.5D+00*(G(2,3)+G(MAXXM,3))
      CALL BOUND(G22)
      CALL XSWEEP(ITER,SUB)
C------ CONVERGENCE RATE CRITERION
      IF(DABS(CIRCO) .LT. 1.0D-05 .OR. DABS(CFST) .LT. 1.0D-06)         -
     1   GO TO 500
      DOG=(CIRC-CIRCO)/CIRCO
      IF(ABS(DOG) .LT. CONVER .AND. DABS(YOR) .LT. CONVER .AND. ITER    -
     1      .GT. 20) GO TO 700
500   CONTINUE
600   CONTINUE
700   CONTINUE
      CALL CPMACH(PP,QQ,NMESH)
      NHNOW=NHNOW+1
C------ SET - UP FOR THE MESH REFINEMENT
      IF(NHNOW .GT. NMESH) GO TO 800
      RLX=RLX*1.04
      CONVER=CONVER*2.
      IF(RLX .GT. 1.95) RLX=1.95
      ITRMAX=TITR2+0.001
      IF(NHNOW .EQ. 3) ITRMAX=TITR3+0.001
      IF(NHNOW .EQ. 4) ITRMAX=TITR*0.5
      TITR=ITRMAX*1.0D+00
      CALL MESH(XCELL,YCELL,MX,NMESH)
      GO TO 100
800   CONTINUE
      RETURN
      END
      SUBROUTINE CONMAP(MAXP,EM,FO,PP,QQ,NTW)
C*************************************************************
C* THIS ROUTINE UTILIZES NEWTON-RAPHSON ITERATIVE PROCEDURE  *
C* TO FIND THE INVERSE FUNCTION Z=Z(W) OF THE BASIC REMAPPING*
C*                        FUNCTION W=W(Z)                    *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      DOUBLE PRECISION HR,VR,PI,TWIST,EM,UEMAX,FMACH,PITCH,P3           -
     2   ,DSIGN,CIRC,GAMMA,AP,AQ,DARCOS,Q4,DVAXY,Q3                     -
     3   ,ZDOL,DF,P4,TO,QO,RO2,RO,BE,BE2,P12,UE,VE,UE1,UE2              -
     4   ,DSQRT,UEC,DVAMXY,ZI,CETR,FO,DVAM,ACRIT
      DOUBLE PRECISION XSI,ETA,EM2,XSI2,ETA2,DVAMX,DVAMY,R1,            -
     2   X2Y2,X2Y2M,D1X,D2X,D3X,D4X,D1,D2,D3,D4,DATAN2,R1Y2,            -
     3   R2,R3,R4,PP,QQ,R1R2,DSIN,DCOS,DLOG,A1,A2,A3,A4,ZR,             -
     4   XSINEW,ETANEW,DIFFX,DIFFY,DABS,XSIO,ETAO,TRI
      DVAM=2.0D+00*EM
      EM2=EM*EM
      DF=1.0D-10
      MPHF=(MAXP+1)/2
      XSI=DCOS(FO)*0.8
      ETA=DSIN(FO)*0.8
      DO 2 I=2,MAXP
      IF(HR(I) .GT. HR(I-1)) GO TO 3
2     CONTINUE
3     CONTINUE
      MID=I-1
      DO 20 K=1,MAXP
      I=K
      IF(NTW .LT. 0) I=MAXP+1-K
      XSIO=HR(I)
      ETAO=VR(I)
      ITER=1
4     CONTINUE
      XSI2=XSI*XSI
      ETA2=ETA*ETA
      DVAMX=DVAM*XSI
      DVAMY=-DVAM*ETA
      DVAXY=-2.0D+00*XSI*ETA
      DVAMXY=DVAXY*EM2
      X2Y2=XSI2+ETA2
      X2Y2M=-XSI2+ETA2
      D1X=EM2-X2Y2
      D2X=1.0D+00-EM2*X2Y2
      D1=DATAN2(DVAMY,D1X)
      IF(D1 .LT. 0.0D+00) D1=D1+2.0D+00*PI
      D2=DATAN2(DVAMY,D2X)
      D3X=-(EM2+X2Y2M)
      D4X=-(1.0D+00+EM2*X2Y2M)
      TRI=D3X*D3X+DVAXY*DVAXY
      CETR=D4X*D4X+DVAMXY*DVAMXY
      R3=DVAM/DSQRT(TRI)
      R4=DVAM/DSQRT(CETR)
      D3=DATAN2(DVAXY,D3X)
      D4=DATAN2(DVAMXY,D4X)
      P3=DSIN(D3)
      Q3=DCOS(D3)
      P4=DSIN(D4)
      Q4=DCOS(D4)
      R1Y2=DVAMY*DVAMY
      R1=D1X*D1X+R1Y2
      R2=D2X*D2X+R1Y2
      R1=DSQRT(R1)/(EM2+X2Y2+DVAMX)
      R2=DSQRT(R2)/(1.0D+00+EM2*X2Y2+DVAMX)
      R1R2=R1*R2
      R1R2=DLOG(R1R2)
      A1=QQ*R1R2+PP*(D2-D1)-XSIO
      R1R2=R1/R2
      R1R2=DLOG(R1R2)
      A2=PP*R1R2+QQ*(D1+D2)-ETAO
      A3=R3*(QQ*Q3-PP*P3)+R4*(QQ*Q4+PP*P4)
      A4=R3*(PP*Q3+QQ*P3)-R4*(PP*Q4-QQ*P4)
      ZDOL=A3*A3+A4*A4
      DIFFX=-(A1*A3+A2*A4)/ZDOL
      DIFFY=-(A2*A3-A1*A4)/ZDOL
      XSI=XSI+DIFFX
      ETA=ETA+DIFFY
      IF(DABS(DIFFX) .LT. DF .AND. DABS(DIFFY) .LT. DF) GO TO 10
      IF(ITER .EQ. 200) GO TO 6
      ITER=ITER+1
      GO TO 4
6     WRITE(6,8)
8     FORMAT(/,20X,27HFAILS TO CONVERGE IN CONMAP,/)
      STOP
10    CONTINUE
      IF(DABS(XSI) .LT. 1.0D-12) GO TO 13
      BE=1.0D+00+X2Y2/EM2
      BE2=BE*BE-4.0D+00*XSI2/EM2
      P12=0.5D+00*(BE-DSQRT(BE2))
      P12=DSQRT(P12)
      IF(P12 .GT. 1.0D+00) P12=1.0D+00
      P12=-DSIGN(P12,XSI)
      VE=DARCOS(P12)-2.0D+00*PI
      IF(ETA .GT. 0.0D+00) VE=-VE-2.0D+00*PI
      IF(ETA .LT. 0.0D+00 .AND. I .GT. MID) VE=VE+2.0D+00*PI
      VE=VE+PI
      UE1=XSI/(EM*P12)
      GO TO 15
13    CONTINUE
      VE=0.5D+00*PI*DSIGN(1.0D+00,VE)
      UE1=1.0D+00+X2Y2/EM2
      UE1=DSQRT(UE1)
15    CONTINUE
      UE2=UE1*UE1-1.0D+00
      UE=DABS(UE1)+DSQRT(UE2)
      S(I)=VE
      F(I)=DLOG(UE)
      IF(I .EQ. 1 .OR. I .EQ. MID .OR. I .EQ. MAXP) GO TO 18
      GO TO 19
18    CONTINUE
      RO=0.95*DSQRT(X2Y2)
      D1=0.05*NTW+DATAN2(ETA,XSI)
      XSI=DCOS(D1)*RO
      ETA=DSIN(D1)*RO
19    CONTINUE
20    CONTINUE
      IF(DABS(PP) .GT. 1.0D-07) GO TO 30
      S(1)=-PI
      S(MAXP)=PI
30    CONTINUE
      F(1)=0.5D+00*(F(1)+F(MAXP))
      F(MAXP)=F(1)
      RETURN
      END
      SUBROUTINE REMAP(GR,ALM,PP,QQ,XC,YC,EM,DRD,NTW,AC,BC,RTE)
C*************************************************************
C* ROUTINE THAT MAPS UNIFORM MESH IN COMPUTATIONAL SPACE BACK*
C*                      INTO THE PHYSICAL SPACE              *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DOUBLE PRECISION UE,DCOSH,VE,DSQRT,DSINH,GR,EM,PP,QQ,DX1,XO,ANA   -
     1 ,AC,BC,DLOG,DATAN2,T1,T2,DSIN,DCOS,UIN,DRD,PI,HR,VR,WV,DA,FMACH  -
     2 ,R1,XC,YC,COEF,US,ALM,EM2,SIS,XSI,ETA,DVAMX,DVAMY,YO,YOY,FL,RTE  -
     3 ,AP,AQ,CIRC,GAMMA,DABS,RR,R2,ETA2,DX2,PITCH,XSI2,SLO,B1,B2,HH
      EM2=EM*EM
      COEF=4.0D+00*PI/XC
      SLO=2.0D+00*PI/PITCH
      B1=ALM/SLO
      B2=PI*QQ/SLO
      NN=1-NTW
      N2N=2*NTW
      DA=2.0D+00*EM
      HH=0.90D+00
      IF(RTE .LE. 0.0001) HH=-0.10-2.*DRD
      HH=DABS(HH)
      IRHS=IMID+1
      ILHS=IMID-1
      IRHSTE=MAXXM-1
      ILHSTE=3
      DO 100 J=2,MAXY
      I1=2
      I2=IMID
      I3=IMID+1
      I4=MAXXM
      IF(J .NE. 2) GO TO 20
      I1=ILHSTE
      I2=ILHS
      I3=IRHS
      I4=IRHSTE
20    CONTINUE
      US=((J-2)*1.0D+00)/YC
      UIN=US*PI
      US=US+(0.38/PI)*DSIN(UIN)
      ANA=US+(HH/PI)*DSIN(UIN)
      UIN=DRD*ANA
      IF(NTW .GT. 0) UIN=-PI+ANA*(PI+DRD)
      DO 60 N=1,2
      DO 50 I=I1,I2
      SIS=COEF*(I-2)
      WV=DSIN(SIS)
      VE=UIN+0.5D+00*SIS+GR*WV
      UE=FI(I)*US
      XSI=EM*DCOSH(UE)*DCOS(VE)
      ETA=EM*DSINH(UE)*DSIN(VE)
      DVAMX=DA*XSI
      DVAMY=-DA*ETA
      XSI2=XSI*XSI
      ETA2=ETA*ETA
      X2Y2=XSI2+ETA2
      DX1=EM2-X2Y2
      DX2=1.0D+00-EM2*X2Y2
      R1=DX1*DX1+DVAMY*DVAMY
      R1=DSQRT(R1)/(EM2+X2Y2+DVAMX)
      R2=DX2*DX2+DVAMY*DVAMY
      R2=DSQRT(R2)/(1.0D+00+EM2*X2Y2+DVAMX)
      T1=DATAN2(DVAMY,DX1)
      IF(T1 .LT. 0.0D+00) T1=2.0D+00*PI+T1
      M=I/IMID
      IF(J .EQ. 2) T1=PI*(NN+N2N*FLOAT(M))
      T2=DATAN2(DVAMY,DX2)
      RR=R1*R2
      RR=DABS(RR)
      XO=QQ*DLOG(RR)+PP*(T2-T1)
      RR=R1/R2
      RR=DABS(RR)
      YO=PP*DLOG(RR)+QQ*(T1+T2)
      XO=XO/SLO+B1
      YO=YO/SLO-B2
      X(I,J)=BC*XO+AC*YO
      Y(I,J)=BC*YO-AC*XO
50    CONTINUE
      I1=I3
      I2=I4
60    CONTINUE
100   CONTINUE
      RETURN
      END
      SUBROUTINE XYINF
C*************************************************************
C*     ROUTINE THAT DEFINES GRID AT AXIAL INFINITIES         *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DOUBLE PRECISION HR,VR,PITCH,POTCH,FMACH,PI,GAMMA,CIRC,AP,AQ
C------ VALUES OF COORDINATES IN DUMMY COLUMNS
      DO 5 J=2,MAXY
      X(1,J)=X(MAXX-2,J)
      X(MAXX,J)=X(3,J)
      Y(1,J)=Y(MAXX-2,J)
      Y(MAXX,J)=Y(3,J)
5     CONTINUE
C------ DELIBERATE MODIFICATIONS - POINTS AT INFINITIES
      IAO=MAXXM-10
      DO 10 I=10,IAO
      IF(X(I,3) .GT. X(I-1,3)) GO TO 12
10    CONTINUE
12    CONTINUE
      IL=I-1
      XMID=X(IL,3)*2.-X(IL,4)
      XILHS=X(ILHS,2)
      XEND=AMIN1(XMID,XILHS)
      X(IMID,2)=XEND
      X(ILHS,2)=XEND
      X(IRHS,2)=XEND
      DO 15 I=IAO,MAXX
      IF(X(I,3) .LT. X(I-1,3)) GO TO 16
15    CONTINUE
16    CONTINUE
      IR=I-1
      X22=X(IR,3)
      DO 20 I=1,10
      II=11-I
      IF(X(II,3) .GT. X22) IR=II
20    CONTINUE
      XMID=X(IR,3)*2.-X(IR,4)
      X32=X(3,2)
      XEND=AMAX1(XMID,X32)
      X(2,2)=XEND
      X(MAXXM,2)=XEND
      X(3,2)=XEND
      X(IRHSTE,2)=XEND
      X(1,2)=XEND
      X(MAXX,2)=XEND
      Y(MAXXM,2)=0.5D+00*(Y(MAXXM,3)+Y(2,3))
      Y(2,2)=Y(MAXXM,2)
      Y(ILHS,2)=2.0D+00*Y(ILHS-1,2)-Y(ILHS-2,2)
      Y(IRHS,2)=Y(ILHS,2)+PITCH
      Y(IMID,2)=0.5D+00*(Y(ILHS,2)+Y(IRHS,2))
C------ DOUBLE VALUED PNTS ON THE CUT
      DO 150 J=2,MAXY
      X2M=0.5D+00*(X(2,J)+X(MAXXM,J))
      X(2,J)=X2M
      X(MAXXM,J)=X2M
      Y2M=0.5D+00*(Y(2,J)+Y(MAXXM,J))
      Y(2,J)=Y2M
      Y(MAXXM,J)=Y2M
150   CONTINUE
      POTCH=0.5D+00*PITCH
      Y(3,2)=Y(2,2)-POTCH
      Y(MAXXM-1,2)=Y(MAXXM,2)+POTCH
C------ VALUES OF COORDINATES IN IMMAGINARY ROWS
      DO 250 I=1,MAXX
      X(I,MAXYP)=2.0D+00*X(I,MAXY)-X(I,MAXYM)
      Y(I,MAXYP)=2.0D+00*Y(I,MAXY)-Y(I,MAXYM)
250   CONTINUE
      DO 280 I=1,IMID
      II=MAXX+1-I
      X(I,1)=X(II,3)
      X(II,1)=X(I,3)
      Y(I,1)=Y(II,3)-PITCH
      Y(II,1)=Y(I,3)+PITCH
280   CONTINUE
      X(2,1)=X(2,2)+PITCH*1.5
      X(IMID,1)=X(IMID,2)-PITCH*1.5
      X(MAXXM,1)=X(2,1)
      X(1,1)=X(3,1)
      X(MAXX,1)=X(1,1)
      Y(IMID,1)=Y(IMID,2)
      Y(2,1)=Y(2,2)
      Y(MAXXM,1)=Y(2,1)
      Y(1,1)=-Y(1,1)
      Y(MAXX,1)=-Y(MAXX,1)
      Y(1,2)=Y(MAXXM-1,2)
      Y(MAXX,2)=Y(3,2)
      RETURN
      END
      SUBROUTINE MESH(XCELL,YCELL,MX,NHALF)
C*************************************************************
C* THIS ROUTINE DOUBLES THE NUMBER OF MESH CELLS IN EACH     *
C* COMPUTATIONAL DIRECTION AND LINEARLY INTERPOLATES VALUES  *
C*       OF THE REDUCED POTENTIAL AT THESE NEW MESH POINTS   *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DOUBLE PRECISION HR,VR,PITCH,POTCH,FMACH,PI,DRD,WH,SIS            -
     2   ,XCELL,YCELL,AP,AQ,CIRC,GAMMA,COEF,DSIN,PP
      I=2
      NON=1
      IF(NHNOW .EQ. 2 .AND. NHALF .EQ. 3) NON=2
      IF(NHNOW .EQ. 2 .AND. NHALF .EQ. 4) NON=4
      IF(NHNOW .EQ. 3 .AND. NHALF .EQ. 4) NON=2
      DO 3 II=2,MX,NON
      FI(I)=SRF(II)
      I=I+1
3     CONTINUE
      XCELL=XCELL*2.0D+00
      YCELL=YCELL*2.0D+00
      MAXXM=XCELL+2.0001
      MAXY=YCELL+2.0001
      MAXX=MAXXM+1
      MAXYP=MAXY+1
      MAXYM=MAXY-1
      IMID=(MAXX+1)/2
      MXXM=MAXXM-2
      MYYM=MAXY-2
C------ RENUMBERING-SHIFTING WITHIN MATRICES
      DO 10 J=2,MAXY,2
      JJ=MAXY+2-J
      JJJ=(JJ+2)/2
      DO 5 I=2,MAXXM,2
      II=MAXXM+2-I
      III=(II+2)/2
      G(II,JJ)=G(III,JJJ)
5     CONTINUE
10    CONTINUE
C------ INTERPOLATION AT MIDPOINTS ALONG 'OLD' HORIZONTALS
      DO 40 J=2,MAXY,2
      DO 30 I=2,MXXM,2
      II=I+1
      G(II,J)=0.5D+00*(G(I,J)+G(I+2,J))
30    CONTINUE
40    CONTINUE
C------ INTERPOLATION AT MIDPOINTS ALONG 'OLD' VERTICALS
      DO 60 I=2,MAXXM,2
      DO 50 J=2,MYYM,2
      JJ=J+1
      G(I,JJ)=0.5D+00*(G(I,J)+G(I,J+2))
50    CONTINUE
60    CONTINUE
C------ INTERPOLATION AT THE CENTER OF THE EACH OLD MESH CELL
      DO 80 I=2,MXXM,2
      II=I+1
      DO 70 J=2,MYYM,2
      JJ=J+1
      G(II,JJ)=0.25D+00*(G(II,J)+G(I,JJ)+G(II,J+2)+G(I+2,JJ))
70    CONTINUE
80    CONTINUE
      DO 100 I=2,IMID
      II=MAXXM+2-I
      G(I,1)=G(II,3)
      G(II,1)=G(I,3)
100   CONTINUE
      G(3,2)=G(2,2)+0.5D+00*CIRC
      G(MAXXM-1,2)=G(MAXXM,2)-0.5D+00*CIRC
      G(IMID+1,2)=0.0D+00
      G(IMID-1,2)=0.0D+00
      RETURN
      END
      SUBROUTINE SPLIF(MUNA,NUNA,KM,VM,KN,VN,MODE,FQM)
C*************************************************************
C*     MODIFIED ROUTINE FROM: SUPERCRITICAL WING SECTIONS II *
C*     CUBIC SPLINE IN ONE SPATIAL COORD. IS FITTED TO DATA  *
C*     ARRAY "F" AT NODES "S" FROM INDEX "M" TO INDEX"N"     *
C*************************************************************
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      DOUBLE PRECISION PI,RADIAN,ALPHA,A2O,CIR,AQ,AP,CIRC,              -
     1YCELL,XCELL,FMACH,SUB,UEMAX,HR,VR,EM,EL,D,C,FPPPP,FF,             -
     2PITCH,GAMMA,B
      IND=0
      M=MUNA
      N=NUNA
      K=IABS(N-M)
      IF (K-1) 81,81,1
1     K=(N-M)/K
      I=M
      J=M+K
      DS=S(J)-S(I)
      D=DS
      IF (DS) 11,81,11
11    DF =(F(J)-F(I))/DS
      IF (KM-2) 12,13,14
12    U=0.5
      V=3.0D+00*(DF-VM)/DS
      GO TO 25
13    U=0.
      V=VM
      GO TO 25
14    U=-1.
      V=-DS*VM
      GO TO 25
21    I=J
      J=J+K
      DS=S(J)-S(I)
      IF (D*DS) 81,81,23
23    DF=(F(J)-F(I))/DS
      B=1.0D+00/(DS+DS+U)
      U=B*DS
      V=B*(6.0D+00*DF-V)
25    FP(I)=U
      FPP(I)=V
      U=(2.-U)*DS
      V=6.0D+00*DF+DS*V
      IF (J-N) 21,31,21
31    IF (KN-2) 32,33,34
32    V=(6.0D+00*VN-V)/U
      GO TO 35
33    V=VN
      GO TO 35
34    V=(DS*VN+FPP(I))/(1.0D+00+FP(I))
35    B=V
41    DS=S(J)-S(I)
      U=FPP(I)-FP(I)*V
      FPPP(I)=(V-U)/DS
      FPP(I)=U
      FP(I)=(F(J)-F(I))/DS-DS*(V+U+U)/6.0D+00
      V=U
      J=I
      I=I-K
      IF(J-M) 41,51,41
51    I=N-K
      FPPP(N)=FPPP(I)
      FPP(N)=SNGL(B)
      FP(N)=DF+D*(FPP(I)+B+B)/6.0D+00
      IND=1
      IF (MODE) 81,81,61
61    FPPP(J)=FQM
      V=FPP(J)
71    I=J
      J=J+K
      DS=S(J)-S(I)
      U=FPP(J)
      FPPP(J)=FPPP(I)+0.5D+00*DS*(F(I)+F(J)-DS*DS*(U+V)/12.0)
      V=U
      IF (J-N) 71,81,71
81    RETURN
      END
      SUBROUTINE INTPL (MI,NI,MUNA,NUNA,MODE)
C*************************************************************
C*     MODIFIED ROUTINE FROM: SUPERCRITICAL WING SECTIONS II *
C*         INTERPOLATION USING PIECEWISE TAYLOR SERIES       *
C*************************************************************
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      DOUBLE PRECISION PI,RADIAN,ALPHA,A2O,CIR,AQ,AP,CIRC,              -
     1YCELL,XCELL,FMACH,SUB,UEMAX,HR,VR,EM,EL,D,C,FPPPP,FF,             -
     2PITCH,GAMMA
      M=MUNA
      N=NUNA
      K=IABS(N-M)
      K=(N-M)/K
      I=M
      MIN=MI
      NIN=NI
      D=S(N)-S(M)
      IF (D*(SI(NI)-SI(MI))) 11,13,13
11    MIN=NI
      NIN=MI
13    KI=IABS(NIN-MIN)
      IF (KI) 21,21,15
15    KI= (NIN-MIN)/KI
21    II=MIN-KI
      C=0.0D+00
      IF (MODE) 31,31,23
23    C=1.0D+00
31    II=II+KI
      SS=SI(II)*1.0D+00
33    I=I+K
      IF (I-N) 35,37,35
35    IF (D*(S(I)-SS)) 33,33,37
37    J=I
      I=I-K
      SS=SS-S(I)
      FPPPP=C*(FPPP(J)-FPPP(I))/(S(J)-S(I))
      FF=FPPP(I)+0.25D+00*SS*FPPPP
      FF=FPP(I)+SS*FF/3.0D+00
      FF=FP(I)+0.5D+00*SS*FF
      FI(II)=F(I)+SS*FF
      IF (II-NIN) 31,41,31
41    RETURN
      END
      SUBROUTINE XSWEEP(ITER,RLX)
C*************************************************************
C*                     FULLY - CONSERVATIVE SCHEME           *
C* SUBROUTINE WHICH SWEEPS COMPUTATIONAL DOMAIN FROM THE L.E.*
C* TO THE RIGHT TO THE T.E. AND THEN AGAIN FROM THE L.E. TO  *
C* THE LEFT TO THE T.E. USING SLOR TECHNIQUE . JAMESON'S     *
C* ROTATED LOCALLY TYPE DEPENDENT FINITE DIFFERENCING WAS    *
C* USED TO EVALUATE THE COEFFICIENTS OF THE CORRECTION MATRIX*
C* RESIDUALS WERE EVALUATED USING FINITE AREA TECHNIQUE AND  *
C* ARTIFICIAL VISCOSITY WAS ADDED IN CONSERVATIVE FORM       *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DIMENSION HRM(2,195),VRM(2,195)
      DOUBLE PRECISION AP,AQ,SUB,CIRC,HR,VR,PI,FMACH,PITCH,A,AI2,       -
     1DSQRT,REOLD,RE,CG,SXX,SYY,SXY,RXX,RYY,RXY,BM,BB,BP,TE,GAMMA,      -
     2DABS,Q2P,REZ,YBETM,YBETP,FYT,RLX
      IR=0
      JR=0
      IG=0
      JG=0
      DOFR=0.
      DOFG=0.
      NS=0
C------ DETERMINATION OF THE LEADING EDGE STAGNATION POINT
      ISTG=ILHS
      IF(FPP(IMID) .LT. FPP(ILHS)) ISTG=IMID
      IF(FPP(IRHS) .LT. FPP(IMID) .AND. ISTG .EQ. IMID) ISTG=IRHS
      IF(ITER .NE. 1) GO TO 5
      WRITE(6,3)
3     FORMAT(45X,25HFULLY-CONSERVATIVE SCHEME,/)
      DO 4 I=1,MAXX
      HRM(1,I)=0.0D+00
      VRM(1,I)=0.0D+00
      HRM(2,I)=0.0D+00
      VRM(2,I)=0.0D+00
      S(I)=0.0D+00
      F(I)=0.0D+00
4     CONTINUE
      A=0.5D+00*(GAMMA-1.0D+00)
      FM2=FMACH*FMACH
      AI2=1.0D+00/FM2+A
      AM=A*FM2
      EX=0.5D+00/A
      ISTG=IMID
      Q2PP=1.
      UIPP=1.
      A2PP=1.
      ROPP=1.
      DETPP=1.
      RXXPP=0.
      RYYPP=0.
      RXYPP=0.
      SXXPP=0.
      SYYPP=0.
      SXYPP=0.
      MAXX2=MAXXM+2
      MAXY2=MAXY+2
5     CONTINUE
      SUB=RLX
      SUBP=RLX
      I1=ISTG-1
      I2=MAXXM
      IST=I1
      IEND=I2
      IMNS=-1
C------ SAVE OLD POTENTIAL FOR THE LEFT SWEEP START
      DO 8 J=1,MAXYP
      HR(J)=0.0D+00
      VR(J)=0.0D+00
      FPPP(J)=0.
      FP(J)=G(ISTG,J)
      GOLD(J)=G(I1-1,J)
8     CONTINUE
10    CONTINUE
C------ X - LOCAL COORDINATE : POINTS LEFTWARD
C------ Y - LOCAL COORDINATE : POINTS DOWNWARD
C------ FIRST SWEEP TO THE RIGHT FROM THE L.E.
      IMNS2=-2*IMNS
      I22=I2+2
      I1P=(1+IMNS)/2
      I1M=(1-IMNS)/2
      DO 1000 II=I1,I2
      I=II*I1M+(I22-II)*I1P
      IP=I-1
      IM=I+1
      IPLUS=I-IMNS
      IMP=I+IMNS
C------ UPSTREAM AND DOWNSTREAM INFINITY ' WINDOWS '
      MO=2
      IF(I .LE. 3 .OR. I .GE. IRHSTE) MO=3
      IF(I .LE. IRHS .AND. I .GE. ILHS) MO=3
      MOM=MO-1
      IJK=MAXX2-I
      HR(1)=HRM(1,IJK)
      VR(1)=VRM(1,IJK)
      YBET=0.
      GIPJPP=G(IPLUS,2)
      GIPJP=G(IPLUS,1)
      GIJPP=G(I,2)
      GIJP=G(I,1)
      XIPJPP=X(IPLUS,2)
      XIPJP=X(IPLUS,1)
      XIJPP=X(I,2)
      XIJP=X(I,1)
      YIPJPP=Y(IPLUS,2)
      YIPJP=Y(IPLUS,1)
      YIJPP=Y(I,2)
      YIJP=Y(I,1)
      XIMJP=X(IM,1)
      XPJP=X(IP,1)
      YIMJP=Y(IM,1)
      YPJP=Y(IP,1)
      MY=MAXY
      IF(I .EQ. IEND .AND. NHNOW .EQ. 1 .AND. ITER .LE. 75) MY=MAXYM
      MYP=MY+1
      HR(2)=HR(2)*(MOM-1)
      VR(2)=VR(2)*(MOM-1)
      XBETP=0.0D+00
      YBETP=0.0D+00
C------ COLLUMN SWEEP BEGINS HERE
      DO 500 J=1,MY
      JP=J+1
      JM=J-1
      JPP=J+2
      JSURF=1-J/MAXY
      JSURFM=1-J/MAXYM
C------ DESCENDING OF ALL PARAMETERS FROM (JP) TO (J)
      SUB=SUBP
      Q2P=Q2PP
      A2P=A2PP
      ROP=ROPP
      DETP=DETPP
      UP=UIPP
      VP=VIPP
      RXX=RXXPP
      RYY=RYYPP
      RXY=RXYPP
      SXX=SXXPP
      SYY=SYYPP
      SXY=SXYPP
      YBETM=YBET
      YBET=YBETP
      XBET=XBETP
      RUX=(RUJP-S(J))
      RV2=RV1
      S(J)=RUJP
      IF(J .EQ. MAXY) GO TO 300
      GIPJ=GIPJP
      GIPJP=GIPJPP
      GIJ=GIJP
      GIJP=GIJPP
      XIJ=XIJP
      XIJP=XIJPP
      XIPJ=XIPJP
      XIPJP=XIPJPP
      YIJ=YIJP
      YIJP=YIJPP
      YIPJ=YIPJP
      YIPJP=YIPJPP
      XIMJ=XIMJP
      XPJ=XPJP
      YIMJ=YIMJP
      YPJ=YPJP
      GIPJPP=G(IPLUS,JPP)
      GIJPP=G(I,JPP)
      XIPJPP=X(IPLUS,JPP)
      XIJPP=X(I,JPP)
      XIMJP=X(IM,JP)
      XPJP=X(IP,JP)
      YIPJPP=Y(IPLUS,JPP)
      YIJPP=Y(I,JPP)
      YIMJP=Y(IM,JP)
      YPJP=Y(IP,JP)
C------ POINT IN THE SWEEP DIRECTION : (I-IMNS/2,JP)
      XX=(XIPJP-XIJP)
      YX=(YIPJP-YIJP)
      XY=XIJPP+XIPJPP-XIJ-XIPJ
      YY=YIJPP+YIPJPP-YIJ-YIPJ
      DET=XX*YY-YX*XY
      GX=(GIPJP-GIJP)
      GY=GIJPP+GIPJPP-GIJ-GIPJ
      UO=AQ+(YY*GX-YX*GY)/DET
      VO=AP+(XX*GY-GX*XY)/DET
      QOJP=UO*UO+VO*VO
      UIJP=(YY*UO-XY*VO)*IMNS
      RUJP=1.0D+00+AM*(1.0D+00-QOJP)
      RUJP=UIJP*(RUJP**EX)
C------ RUJP SHOULD BE DEVIDED BY 8.
C------ LOWER CELL : (J+1/2)
      XX=XPJ+XPJP-XIMJ-XIMJP
      YX=YPJ+YPJP-YIMJ-YIMJP
      XY=XIJP-XIJ
      YY=YIJP-YIJ
      DET=XX*YY-YX*XY
      GX=(GIPJ+GIPJP-GOLD(J)-GOLD(JP))*IMNS
      GY=GIJP-GIJ
      UO=AQ+(YY*GX-YX*GY)/DET
      VO=AP+(XX*GY-GX*XY)/DET
      QO1=UO*UO+VO*VO
      VI1=(XX*VO-YX*UO)
      RV1=1.0D+00+AM*(1.0D+00-QO1)
      RV1=VI1*(RV1**EX)
C------ RV1 SHOULD BE DEVIDED BY 8.
C------ CENTRAL CELL - POINT JP=J+1
      XX=0.25D+00*(XPJP-XIMJP)
      YX=0.25D+00*(YPJP-YIMJP)
      XY=0.25D+00*(XIJPP-XIJ)
      YY=0.25D+00*(YIJPP-YIJ)
      DETPP=XX*YY-XY*YX
      XX=XX/DETPP
      YX=YX/DETPP
      XY=XY/DETPP
      YY=YY/DETPP
      GX=0.25D+00*(GIPJP-GOLD(JP))*IMNS
      GY=0.25D+00*(GIJPP-GIJ)
C------ PHYSICAL VELOCITY COMPONENTS
      UO=AQ+(YY*GX-YX*GY)
      VO=AP+(XX*GY-GX*XY)
C------ MODIFIED CONTRAVARIANT VELOCITY COMPONENTS
      UIPP=YY*UO-XY*VO
      VIPP=(XX*VO-YX*UO)*JSURFM
      Q2PP=UO*UO+VO*VO
      A2PP=AI2-A*Q2PP
      ROPP=A2PP*FM2
      ROPP=ROPP**EX
      SXXPP=ROPP*((XY*XY+YY*YY)-UIPP*UIPP/A2PP)
      SYYPP=ROPP*((XX*XX+YX*YX)-VIPP*VIPP/A2PP)
      SXYPP=ROPP*((XX*XY+YX*YY)+UIPP*VIPP/A2PP)*IMNS2
      XBETP=0.0D+00
      YBETP=0.0D+00
      SUBP=RLX
      RXXPP=0.0D+00
      RYYPP=0.0D+00
      RXYPP=0.0D+00
C------ ARTIFICIAL VISCOSITY TERMS
      IF(Q2PP .LT. A2PP) GO TO 300
      NS=NS+1
      SUBP=2.0D+00
      SWPP=ROPP*(1.0D+00/A2PP-1.0D+00/Q2PP)
      FXX=GIPJP-2.0D+00*GIJP+GOLD(JP)
      FYY=GIJPP-2.0D+00*GIJP+GIJ
      FXY=(GIPJPP+GOLD(J)-GOLD(JPP)-GIPJ)*0.125D+00
      RXXPP=-SWPP*UIPP*UIPP
      RYYPP=-SWPP*VIPP*VIPP
      RXYPP=IMNS2*SWPP*UIPP*VIPP
      XBETP=(RXXPP*FXX+RXYPP*FXY)*DETPP
      YBETP=(RYYPP*FYY+RXYPP*FXY)*DETPP
300   CONTINUE
      IF(J .LE. MOM .OR. I .EQ. IST) GO TO 500
C------ ELEMENTS OF THE CORRECTION MATRIX
      JMNS=1
      IF(VP .LT. (-1.0D-08)) JMNS=-1
      JJP=(1+JMNS)/2
      JJM=(1-JMNS)/2
      RXYIJ=JMNS*RXY
C------ ELEMENTS OF THE FINAL CORRECTION MATRIX
      BM=-SYY+RYY*(2+JMNS)+JJP*RXY
      BB=2.0D+00*(SYY-RXYIJ+SXX/SUB)-3.0D+00*(RXX+RYY)-RYY
      BM=-SYY+RYY*(2-JMNS)-JJM*RXY
C------ ARTIFICIAL VISCOSITY
      RVISC=(JJM*(YBET-YBETP)+JJP*(YBET-YBETM))
      RVISC=(RVISC+(XBET-FPPP(J)))/DETP
      FPPP(J)=XBET
C------ TOTAL RESIDUE
      RE=(0.25D+00/DETP)*(RUX+RV1*JSURF+(JSURF-2)*RV2)
      REOLD=(3.0D+00*RXX-SXX+RXYIJ)*(GOLD(J)-G(IMP,J))
      FYT=GOLD(JP)-G(IMP,JP)-GOLD(JM)+G(IMP,JM)
      REOLD=REOLD+0.25D+00*(SXY-RXY)*FYT
      REZ=RE+REOLD-RVISC
      IF(DABS(REZ) .LT. ABS(DOFR)) GO TO 450
      DOFR=REZ
      IR=I
      JR=J
450   CONTINUE
C------ INVERSION OF THE RELAXATION MATRIX
      BD=1.0D+00/(BB-BM*HR(JM))
      HR(J)=BD*BP
      VR(J)=BD*(REZ-BM*VR(JM))
500   CONTINUE
      HRM(2,IJK)=HR(3)
      VRM(2,IJK)=VR(3)
      FPP(I)=Q2P
      IF(I .NE. IST) GO TO 540
      DO 525 J=1,MAXYP
      SI(J)=-S(J)
      FI(J)=F(J)
525   CONTINUE
540   CONTINUE
C------ CORRECTION TO THE POTENTIAL
C------ SAVE OLD VALUES OF THE POTENTIAL
      CG=0.0D+00
      GOLD(1)=G(I,1)
      GOLD(MOM)=G(I,MOM)
      GOLD(MYP)=G(I,MYP)
      J1=2+(MAXY-MY)
      J2=MAXYP-MOM
      DO 600 JJ=J1,J2
      J=MAXY2-JJ
      CG=-VR(J)-HR(J)*CG
      IF(DABS(CG) .LE. ABS(DOFG)) GO TO 550
      DOFG=CG
      IG=I
      JG=J
550   CONTINUE
      GOLD(J)=G(I,J)
      G(I,J)=GOLD(J)-CG
600   CONTINUE
1000  CONTINUE
      IF(I .EQ. 2) GO TO 1010
C------ GET READY FOR THE LEFT SWEEP
      DO 1005 J=1,MAXYP
      FPPP(J)=0.0D+00
      GOLD(J)=FP(J)
      S(J)=SI(J)
      F(J)=FI(J)
1005  CONTINUE
      HR(2)=0.0D+00
      VR(2)=0.0D+00
      I2=I1
      I1=2
      IEND=2
      IST=1
      IMNS=1
      GO TO 10
1010  CONTINUE
      DO 1025 I=2,MAXXM
      HRM(1,I)=HRM(2,I)
      VRM(1,I)=VRM(2,I)
1025  CONTINUE
      IF( MY .EQ. MAXY) GO TO 1050
      G(MAXXM,MAXY)=0.5D+00*(G(IRHSTE,MAXY)+G(MAXXM,MAXYM))
      G(2,MAXY)=0.5D+00*(G(3,MAXY)+G(2,MAXYM))
 1050 CONTINUE
C------ CIRCULATION
      CIRC=G(MAXXM,MAXY)-G(2,MAXY)
      WRITE(6,1100) ITER,IR,JR,DOFR,IG,JG,DOFG,CIRC,SUB,ISTG,NS
1100  FORMAT(8X,2I6,I4,1X,E14.6,I6,I5,1X,3E14.6,I4,I6)
      RETURN
      END
      SUBROUTINE BOUND(G22)
C*************************************************************
C*     ENFORCEMENT OF BOUNDARY AND PERIODICITY CONDITIONS    *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DOUBLE PRECISION HR,VR,PI,FMACH,PITCH,GAMMA,AP,AQ,CIRC,G22,AR
C------ DOWNSTREAM INFINITY B.C.
      G(3,2)=G22
      G(IRHSTE,2)=G22
      G(2,2)=G22-0.5D+00*CIRC
      G(MAXXM,2)=G22+0.5D+00*CIRC
C------ VERTICAL WALLS(POTENTIAL JUMP ACCROSS THE CUT)
      DO 5 J=2,MAXY
      G(1,J)=G(MAXX-2,J)-CIRC
      G(MAXX,J)=G(3,J)+CIRC
      G(2,J)=0.5D+00*(G(2,J)+G(MAXXM,J)-CIRC)
      G(MAXXM,J)=G(2,J)+CIRC
5     CONTINUE
C------ UPSTREAM INFINITY
      G(ILHS,2)=0.0D+00
      G(IMID,2)=0.0D+00
      G(IRHS,2)=0.0D+00
C------ (UPPER WALL - PERIODIC B.C. )
      DO 10 I=3,ILHS
      II=MAXXM+2-I
      G(I,2)=0.5D+00*(G(I,2)+G(II,2))
      G(II,2)=G(I,2)
      G(I,1)=G(II,3)
      G(II,1)=G(I,3)
10    CONTINUE
C------ (LOWER WALL - AIRFOIL SURFACE B.C. )
      DO 20 I=1,MAXX
      G(I,MAXYP)=2.0D+00*G(I,MAXY)-G(I,MAXYM)
20    CONTINUE
      RETURN
      END
      SUBROUTINE CPMACH(PP,QQ,NMESH)
C***********************************************************
C* ROUTINE WHICH GENERATES THE MACH NUMBER CHART,CALCULATES*
C* RELATIVE SURFACE VELOCITY,DENSITY,PRESSURE COEFFICIENTS *
C*                 AND GENERATES THE CP-CHART              *
C***********************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DIMENSION NPT(195)
      DATA BLANK,DOT,STAR,PLUS/1H ,1H.,1H*,1H+/
      DOUBLE PRECISION HR,VR,PI,FMACH,PITCH,CHRDA,DOL,CPC,CPO,          -
     1PP,QQ,DS,SRFM,DSQRT,DABS,AP,AQ,CIRC,GAMMA
      AAA=0.5D+00*(GAMMA-1.0D+00)
      A1A=GAMMA/(GAMMA-1.0D+00)
      DS=1.0D+00/(GAMMA-1.0D+00)
      AI2=1.0D+00/(FMACH*FMACH)+AAA
      DOL=2.0D+00/(GAMMA*FMACH*FMACH)
      CPC=AAA*FMACH*FMACH
      CPCRIT=2.0D+00*(1.0D+00+CPC)/(GAMMA+1.0D+00)
      CPCRIT=DOL*(CPCRIT**A1A-1.0D+00)
      GAM=GAMMA*1.0D+00
      IO=(IMID+MAXXM)/2
      IF(NHNOW .LT. NMESH .OR. FMACH .LT. 0.0011) GO TO 41
      WRITE(6,5)
5     FORMAT(//,40X,40HMACH NUMBER CHART IN COMPUTATIONAL PLANE,/)
      DO 8 J=2,50
      NPT(J-1)=J
8     CONTINUE
      WRITE(6,10) (NPT(J),J=1,49)
10    FORMAT(6X,49I2,/)
      DO 12 J=2,50
      NPT(J)=0
12    CONTINUE
15    FORMAT(1X,I4,1X,49I2)
      I1=2
      IF(DABS(CIRC) .LT. 1.0D-04) I1=IMID
      DO 25 I=I1,MAXXM
      IP=I+1
      IM=I-1
      DO 20 J=2,MAXY
      JP=J+1
      JM=J-1
      XX=X(IP,J)-X(IM,J)
      YX=Y(IP,J)-Y(IM,J)
      XY=X(I,JP)-X(I,JM)
      YY=Y(I,JP)-Y(I,JM)
      DET=1.0D+00/(XX*YY-XY*YX)
      GX=G(IP,J)-G(IM,J)
      GY=G(I,JP)-G(I,JM)
      UO=AQ+(YY*GX-YX*GY)*DET
      VO=AP+(XX*GY-GX*XY)*DET
      Q2=UO*UO+VO*VO
      A2=AI2-AAA*Q2
      FMCH=Q2/A2
      FMCH=SQRT(FMCH)
      HR(J)=FMCH
      NPT(J)=FMCH*10.
      IF(FMCH .GT. 9.99) NPT(J)=99
      VR(J)=VO/UO
20    CONTINUE
      NPT(1)=I
      WRITE(6,15) (NPT(J),J=1,50)
25    CONTINUE
      ILE=I-1
      EXTM2=0.5D+00*(HR(3)+HR(4))
      WRITE(6,30) EXTM2
30    FORMAT(/,40X,16HEXIT MACH NO. = ,E13.5)
      Q2=AI2*EXTM2*EXTM2/(1.0D+00+AAA*EXTM2*EXTM2)
      RO=1.0D+00+CPC*(1.0D+00-Q2)
      RO=RO**DS
      WRITE(6,33) RO
33    FORMAT(40X,16HEXIT DENSITY  = ,E13.5)
      EXT=0.5D+00*(VR(3)+VR(4))
      EXT=ATAN(EXT)
      EXTANG=EXT*180./PI
      WRITE(6,35) EXTANG
35    FORMAT(40X,16HEXIT ANGLE    = ,E13.5)
C------ CRITICAL COEFFICIENT OF PRESSURE
      CP=2.0D+00*(1.0D+00+CPC)/(GAMMA+1.0D+00)
      CP=DOL*(CP**A1A-1.0D+00)
      WRITE(6,36) CP
36    FORMAT(40X,16HC.P. CRITICAL = ,E13.5)
      CPO=DOL*((1.0D+00+CPC)**A1A-1.0D+00)
      WRITE(6,40) CPO
40    FORMAT(40X,16HCP-STAGNATION = ,E13.5)
41    CONTINUE
      WRITE(6,42)
42    FORMAT(/,38X,31HAIRFOIL SURFACE FLOW PARAMETERS)
      WRITE(6,44)
44    FORMAT(38X,31H-------------------------------)
      WRITE(6,45)
45    FORMAT(/,3X,1HI,7X,1HX,11X,1HY,13X,5HXNORM,13X,2HCP,              -
     1 10X,4HDENS,11X,4HMACH,11X,6HQ/QINF,/)
50    FORMAT(1X,I3,7E14.6)
C------ VELOCITY COMPONENTS ON THE AIRFOIL SURFACE
      DO 60 I=2,MAXXM
      VEL=FPP(I)
      SRFM=VEL/(AI2-AAA*VEL)
C------ LOCAL MACH NUMBER ON THE SURFACE OF THE AIRFOIL
      SRFM=DSQRT(SRFM)
C------ COEFFICIENT OF PRESSURE
      CPO=1.0D+00+CPC*(1.0D+00-VEL)
      DENS=CPO**DS
      CP=DOL*(CPO**A1A-1.0D+00)
      F(I)=CP
      VEL=SQRT(VEL)
      XNORM=-QQ*X(I,MAXY)+PP*Y(I,MAXY)
      XNORM=ABS(XNORM)
      S(I)=XNORM
      WRITE(6,50)  I,X(I,MAXY),Y(I,MAXY),XNORM,CP,DENS,SRFM,VEL
60    CONTINUE
      IF(NHNOW .LT. NMESH) GO TO 101
      WRITE(6,62)
62    FORMAT(/,4X,3HX/C,40X,15HCP-DISTRIBUTION,/)
      CPMIN=F(2)
      CPMAX=F(2)
      DO 65 I=3,MAXXM
      IF(F(I) .GT. CPMIN) CPMIN=F(I)
      IF(F(I) .LT. CPMAX) CPMAX=F(I)
65    CONTINUE
70    FORMAT(F8.4,2X,100A1)
      SCALE=100./(CPMIN-CPMAX)
      NON=CPMIN*SCALE
      DO 75 I=11,110
      FPPP(I)=DOT
75    CONTINUE
      NX=SCALE+0.0001
      NCRT=(CPMIN-CPCRIT)*SCALE
      IF(NCRT .GT. 115) NCRT=1
      DO 100 L=1,NX
      LM=L-1
      FPPP(NCRT)=PLUS
      DO 90 II=2,MAXXM
      N=S(II)*SCALE
      M=1
      IF(N .EQ. LM) M=(CPMIN-F(II))*SCALE
      FPPP(M)=STAR
90    CONTINUE
      FPPP(NON)=DOT
      XNORM=L*(CPMIN-CPMAX)
      WRITE(6,70) XNORM,(FPPP(I),I=11,110)
      DO 95 I=11,110
      FPPP(I)=BLANK
95    CONTINUE
100   CONTINUE
101   CONTINUE
      RETURN
      END
      SUBROUTINE XSWEEP(ITER,RLX)
C*************************************************************
C*                     QUASI - CONSERVATIVE SCHEME           *
C* SUBROUTINE WHICH SWEEPS COMPUTATIONAL DOMAIN FROM THE L.E.*
C* TO THE RIGHT TO THE T.E. AND THEN AGAIN FROM THE L.E. TO  *
C* THE LEFT TO THE T.E. USING SLOR TECHNIQUE . JAMESON'S     *
C* ROTATED LOCALLY TYPE DEPENDENT FINITE DIFFERENCING WAS    *
C* USED TO EVALUATE THE COEFFICIENTS OF THE CORRECTION MATRIX*
C* RESIDUALS WERE EVALUATED USING FINITE AREA TECHNIQUE AND  *
C* ARTIFICIAL VISCOSITY WAS ADDED IN CONSERVATIVE FORM       *
C*************************************************************
      COMMON/BLK1/ HR(195),VR(195),PI,FMACH,PITCH,GAMMA,                -
     1              AP,AQ,CIRC,MAXX,MAXXM,MAXY,MAXYM,MAXYP,             -
     2             IMID,NHNOW,ILHS,IRHS,ILHSTE,IRHSTE
      COMMON/BLK2/ S(195),F(195),SI(195),FI(195),FP(195),FPP(195),      -
     1             FPPP(195),GOLD(195),SRF(195)
      COMMON/BLK3/ X(195,51),Y(195,51),G(195,51)
      DIMENSION HRM(2,195),VRM(2,195)
      DOUBLE PRECISION AP,AQ,SUB,CIRC,HR,VR,PI,FMACH,PITCH,A,AI2,       -
     1DSQRT,REOLD,RE,CG,SXX,SYY,SXY,RXX,RYY,RXY,BM,BB,BP,TE,GAMMA,      -
     2DABS,Q2P,REZ,YBETM,YBETP,FYT,RLX
      IR=0
      JR=0
      IG=0
      JG=0
      DOFR=0.
      DOFG=0.
      NS=0
C------ DETERMINATION OF THE LEADING EDGE STAGNATION POINT
      ISTG=ILHS
      IF(FPP(IMID) .LT. FPP(ILHS)) ISTG=IMID
      IF(FPP(IRHS) .LT. FPP(IMID) .AND. ISTG .EQ. IMID) ISTG=IRHS
      IF(ITER .NE. 1) GO TO 5
      WRITE(6,3)
3     FORMAT(45X,25HQUASI-CONSERVATIVE SCHEME,/)
      DO 4 I=1,MAXX
      HRM(1,I)=0.0D+00
      VRM(1,I)=0.0D+00
      HRM(2,I)=0.0D+00
      VRM(2,I)=0.0D+00
      S(I)=0.0D+00
      F(I)=0.0D+00
4     CONTINUE
      A=0.5D+00*(GAMMA-1.0D+00)
      FM2=FMACH*FMACH
      AI2=1.0D+00/FM2+A
      EX=0.5D+00/A
      ISTG=IMID
      QO1=0.
      VI1=0.
      Q2PP=1.
      UIPP=1.
      VIPP=0.
      A2PP=1.
      ROPP=1.
      DETPP=1.
      YBETP=0.
      XBETP=0.
      RXXPP=0.
      RYYPP=0.
      RXYPP=0.
      SXXPP=0.
      SYYPP=0.
      SXYPP=0.
      UIJP=1.
      QOJP=1.
      MAXX2=MAXXM+2
      MAXY2=MAXY+2
5     CONTINUE
      SUB=RLX
      SUBP=RLX
      I1=ISTG-1
      I2=MAXXM
      IST=I1
      IEND=I2
      IMNS=-1
C------ SAVE OLD POTENTIAL FOR THE LEFT SWEEP START
      DO 8 J=1,MAXYP
      HR(J)=0.0D+00
      VR(J)=0.0D+00
      FPPP(J)=0.
      FP(J)=G(ISTG,J)
      GOLD(J)=G(I1-1,J)
8     CONTINUE
10    CONTINUE
C------ X - LOCAL COORDINATE : POINTS LEFTWARD
C------ Y - LOCAL COORDINATE : POINTS DOWNWARD
C------ FIRST SWEEP TO THE RIGHT FROM THE L.E.
      IMNS2=-2*IMNS
      I22=I2+2
      DO 1000 II=I1,I2
      I=II*1
      IF(IMNS .GT. 0) I=I22-II
      IP=I-1
      IM=I+1
      IPLUS=I-IMNS
      IMP=I+IMNS
C------ UPSTREAM AND DOWNSTREAM INFINITY ' WINDOWS '
      MO=2
      IF(I .LE. 3 .OR. I .GE. IRHSTE) MO=3
      IF(I .LE. IRHS .AND. I .GE. ILHS) MO=3
      MOM=MO-1
      IJK=MAXX2-I
      HR(1)=HRM(1,IJK)
      VR(1)=VRM(1,IJK)
      GIPJPP=G(IPLUS,2)
      GIPJP=G(IPLUS,1)
      GIJPP=G(I,2)
      GIJP=G(I,1)
      XIPJPP=X(IPLUS,2)
      XIPJP=X(IPLUS,1)
      XIJPP=X(I,2)
      XIJP=X(I,1)
      YIPJPP=Y(IPLUS,2)
      YIPJP=Y(IPLUS,1)
      YIJPP=Y(I,2)
      YIJP=Y(I,1)
      XIMJP=X(IM,1)
      XPJP=X(IP,1)
      YIMJP=Y(IM,1)
      YPJP=Y(IP,1)
      MY=MAXY
      IF(I .EQ. IEND .AND. NHNOW .EQ. 1 .AND. ITER .LE. 75) MY=MAXYM
      MYP=MY+1
      HR(2)=HR(2)*(MOM-1)
      VR(2)=VR(2)*(MOM-1)
      XBETP=0.0D+00
      YBETP=0.0D+00
      YBET=0.0D+00
C------ COLLUMN SWEEP BEGINS HERE
      DO 500 J=1,MY
      JP=J+1
      JM=J-1
      JPP=J+2
      JSURF=1-J/MAXY
      JSURFM=1-J/MAXYM
C------ DESCENDING OF ALL PARAMETERS FROM (JP) TO (J)
      SUB=SUBP
      VI2=VI1
      QO2=QO1
      UP=UIPP
      VP=VIPP
      Q2P=Q2PP
      A2P=A2PP
      ROP=ROPP
      DETP=DETPP
      RXX=RXXPP
      RYY=RYYPP
      RXY=RXYPP
      SXX=SXXPP
      SYY=SYYPP
      SXY=SXYPP
      YBETM=YBET
      YBET=YBETP
      XBET=XBETP
      UX=(UIJP-S(J))
      QX=(QOJP-F(J))*IMNS
      S(J)=UIJP
      F(J)=QOJP
      IF(J .EQ. MAXY) GO TO 300
      GIPJ=GIPJP
      GIPJP=GIPJPP
      GIJ=GIJP
      GIJP=GIJPP
      XIJ=XIJP
      XIJP=XIJPP
      XIPJ=XIPJP
      XIPJP=XIPJPP
      YIJ=YIJP
      YIJP=YIJPP
      YIPJ=YIPJP
      YIPJP=YIPJPP
      XIMJ=XIMJP
      XPJ=XPJP
      YIMJ=YIMJP
      YPJ=YPJP
      GIPJPP=G(IPLUS,JPP)
      GIJPP=G(I,JPP)
      XIPJPP=X(IPLUS,JPP)
      XIJPP=X(I,JPP)
      XIMJP=X(IM,JP)
      XPJP=X(IP,JP)
      YIPJPP=Y(IPLUS,JPP)
      YIJPP=Y(I,JPP)
      YIMJP=Y(IM,JP)
      YPJP=Y(IP,JP)
C------ POINT IN THE SWEEP DIRECTION : (I-IMNS/2,JP)
      XX=(XIPJP-XIJP)
      YX=(YIPJP-YIJP)
      XY=XIJPP+XIPJPP-XIJ-XIPJ
      YY=YIJPP+YIPJPP-YIJ-YIPJ
      DET=XX*YY-YX*XY
      GX=(GIPJP-GIJP)
      GY=GIJPP+GIPJPP-GIJ-GIPJ
      UO=AQ+(YY*GX-YX*GY)/DET
      VO=AP+(XX*GY-GX*XY)/DET
      QOJP=UO*UO+VO*VO
      UIJP=(YY*UO-XY*VO)*IMNS
C------ LOWER CELL : (I,J+1/2)
      XX=XPJ+XPJP-XIMJ-XIMJP
      YX=YPJ+YPJP-YIMJ-YIMJP
      XY=XIJP-XIJ
      YY=YIJP-YIJ
      DET=XX*YY-YX*XY
      GX=(GIPJ+GIPJP-GOLD(J)-GOLD(JP))*IMNS
      GY=GIJP-GIJ
      UO=AQ+(YY*GX-YX*GY)/DET
      VO=AP+(XX*GY-GX*XY)/DET
      QO1=UO*UO+VO*VO
      VI1=XX*VO-YX*UO
C------ CENTRAL CELL : (I,J+1)
      XX=0.25D+00*(XPJP-XIMJP)
      YX=0.25D+00*(YPJP-YIMJP)
      XY=0.25D+00*(XIJPP-XIJ)
      YY=0.25D+00*(YIJPP-YIJ)
      DETPP=XX*YY-XY*YX
      XX=XX/DETPP
      YX=YX/DETPP
      XY=XY/DETPP
      YY=YY/DETPP
      GX=0.25D+00*(GIPJP-GOLD(JP))*IMNS
      GY=0.25D+00*(GIJPP-GIJ)
C------ PHYSICAL VELOCITY COMPONENTS
      UO=AQ+(YY*GX-YX*GY)
      VO=AP+(XX*GY-GX*XY)
C------ MODIFIED CONTRAVARIANT VELOCITY COMPONENTS
      UIPP=YY*UO-XY*VO
      VIPP=(XX*VO-YX*UO)*JSURFM
      Q2PP=UO*UO+VO*VO
      A2PP=AI2-A*Q2PP
      ROPP=A2PP*FM2
      ROPP=ROPP**EX
      SXXPP=ROPP*((XY*XY+YY*YY)-UIPP*UIPP/A2PP)
      SYYPP=ROPP*((XX*XX+YX*YX)-VIPP*VIPP/A2PP)
      SXYPP=ROPP*((XX*XY+YX*YY)+UIPP*VIPP/A2PP)*IMNS2
      XBETP=0.0D+00
      YBETP=0.0D+00
      SUBP=RLX
      RXXPP=0.0D+00
      RYYPP=0.0D+00
      RXYPP=0.0D+00
C------ ARTIFICIAL VISCOSITY TERMS
      IF(Q2PP .LT. A2PP) GO TO 300
      NS=NS+1
      SUBP=2.0D+00
      SWPP=ROPP*(1.0D+00/A2PP-1.0D+00/Q2PP)
      FXX=GIPJP-2.0D+00*GIJP+GOLD(JP)
      FYY=GIJPP-2.0D+00*GIJP+GIJ
      FXY=(GIPJPP+GOLD(J)-GOLD(JPP)-GIPJ)*0.125D+00
      RXXPP=-SWPP*UIPP*UIPP
      RYYPP=-SWPP*VIPP*VIPP
      RXYPP=-IMNS2*SWPP*UIPP*VIPP
      XBETP=(RXXPP*FXX+RXYPP*FXY)*DETPP
      YBETP=(RXYPP*FXY+RYYPP*FYY)*DETPP
300   CONTINUE
      IF(J .LE. MOM .OR. I .EQ. IST) GO TO 500
C------ ELEMENTS OF THE CORRECTION MATRIX
      JMNS=1
      IF(VP .LT. (-1.0D-08)) JMNS=-1
      JJP=(1+JMNS)/2
      JJM=(1-JMNS)/2
      RXYIJ=JMNS*RXY
C------ ELEMENTS OF THE FINAL CORRECTION MATRIX
      BM=-SYY+RYY*(2+JMNS)+JJP*RXY
      BB=2.0D+00*(SYY-RXYIJ+SXX/SUB)-3.0D+00*(RXX+RYY)-RYY
      BM=-SYY+RYY*(2-JMNS)-JJM*RXY
C------ ARTIFICIAL VISCOSITY
      RVISC=(JJM*(YBET-YBETP)+JJP*(YBET-YBETM))
      RVISC=(RVISC+(XBET-FPPP(J)))/DETP
      FPPP(J)=XBET
C------ TOTAL RESIDUE
      RE=(0.25D+00/DETP)*(UX+VI1*JSURF+(JSURF-2)*VI2)
      RE=RE-(UP*QX+VP*(QO1-QO2))/A2P
      REOLD=(3.0D+00*RXX-SXX+RXYIJ)*(GOLD(J)-G(IMP,J))
      FYT=GOLD(JP)-G(IMP,JP)-GOLD(JM)+G(IMP,JM)
      REOLD=REOLD+0.25D+00*(SXY-RXY)*FYT
      REZ=ROP*RE+REOLD-RVISC
      IF(DABS(REZ) .LT. ABS(DOFR)) GO TO 450
      DOFR=REZ
      IR=I
      JR=J
450   CONTINUE
C------ INVERSION OF THE RELAXATION MATRIX
      BB=1.0D+00/(BB-BM*HR(JM))
      HR(J)=BB*BP
      VR(J)=BB*(REZ-BM*VR(JM))
500   CONTINUE
      HRM(2,IJK)=HR(3)
      VRM(2,IJK)=VR(3)
      FPP(I)=Q2P
      IF(I .NE. IST) GO TO 540
      DO 525 J=1,MAXYP
      SI(J)=-S(J)
      FI(J)=F(J)
525   CONTINUE
540   CONTINUE
C------ CORRECTION TO THE POTENTIAL
C------ SAVE OLD VALUES OF THE POTENTIAL
      CG=0.0D+00
      GOLD(1)=G(I,1)
      GOLD(MOM)=G(I,MOM)
      GOLD(MYP)=G(I,MYP)
      J1=2+(MAXY-MY)
      J2=MAXYP-MOM
      DO 600 JJ=J1,J2
      J=MAXY2-JJ
      CG=-VR(J)-HR(J)*CG
      IF(DABS(CG) .LE. ABS(DOFG)) GO TO 550
      DOFG=CG
      IG=I
      JG=J
550   CONTINUE
      GOLD(J)=G(I,J)
      G(I,J)=GOLD(J)-CG
600   CONTINUE
1000  CONTINUE
      IF(I .EQ. 2) GO TO 1010
C------ GET READY FOR THE LEFT SWEEP
      DO 1005 J=1,MAXYP
      FPPP(J)=0.0D+00
      GOLD(J)=FP(J)
      S(J)=SI(J)
      F(J)=FI(J)
1005  CONTINUE
      I2=I1
      I1=2
      IEND=2
      IST=1
      IMNS=1
      GO TO 10
1010  CONTINUE
      DO 1025 I=2,MAXXM
      HRM(1,I)=HRM(2,I)
      VRM(1,I)=VRM(2,I)
1025  CONTINUE
      IF( MY .EQ. MAXY) GO TO 1050
      G(MAXXM,MAXY)=0.5D+00*(G(IRHSTE,MAXY)+G(MAXXM,MAXYM))
      G(2,MAXY)=0.5D+00*(G(3,MAXY)+G(2,MAXYM))
 1050 CONTINUE
C------ CIRCULATION
      CIRC=G(MAXXM,MAXY)-G(2,MAXY)
      WRITE(6,1100) ITER,IR,JR,DOFR,IG,JG,DOFG,CIRC,SUB,ISTG,NS
1100  FORMAT(8X,2I6,I4,1X,E14.6,I6,I5,1X,3E14.6,I4,I6)
      RETURN
      END
